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Expansion of natural gas networks is a critical process involving substantial capital expenditures with 
complex decision-support requirements. Given the non-convex nature of gas transmission constraints, global 
optimality and infeasibility guarantees can only be offered by global optimisation approaches. Unfortunately, 
state-of-the-art global optimisation solvers are unable to scale up to real-world size instances. In this study, 
we present a convex mixed-integer second-order cone relaxation for the gas expansion planning problem under 
steady-state conditions. The underlying model offers tight lower bounds with high computational efficiency. 
In addition, the optimal solution of the relaxation can often be used to derive high-quality solutions to the 
original problem, leading to provably tight optimality gaps and, in some cases, global optimal solutions. The 
convex relaxation is based on a few key ideas, including the introduction of flux direction variables, exact 
McCormick relaxations, on/off constraints, and integer cuts. Numerical experiments are conducted on the 
traditional Belgian gas network, as well as other real larger networks. The results demonstrate both the 
accuracy and computational speed of the relaxation and its ability to produce high-quality solutions. 


1. Introduction 

In recent years, the construction of natural gas pipelines has witnessed a tremendous 
growth on a world-wide level. In the U.S., for instance, a $3 billion expansion project of 
the gas pipeline system in New England is planned for late 2016. In Europe, the European 
Investment Bank is supporting a €98 million project for the expansion of gas pipelines in 
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western Poland, to be completed by 2017. These expansion projects aim at increasing gas 
flow capacity on existing pipeline systems and/or bringing new gas wells into production 
and commercialization. In addition, the expansion or reinforcement of a pipeline network 
can also be considered as a risk-awareness strategy to fulfill short or long-term operational 
management requirements when unforeseen events occur such as component failures or 
excessive stress and congestion due to extreme weather conditions. These events were 
observed in New England during the polar vortex experienced in January 2014, when 
major gas-fired power plants in the northeast of the U.S. were forced to shut down due 
to mechanical problems and shortages of gas fuel supplies, which drove wholesale power 
prices up 

According to the U.S. Energy Information Administration (EIA), a project 
for the development and expansion of a Gas Transmission Network (GTN) 
takes an average of three years from its first announcement until its comple¬ 


tion (U.S. Energy Information Administration 2008). The project starts by determining 
the market needs within an open season exercise where nonbinding agreements of capac¬ 
ity rights are offered to potential customers. The second step consists in developing the 
expansion design with initial financial commitments from the potential customers. Note 
that expansions of the gas system include the installation parallel pipelines along existing 
ones (looping), the conversion of oil pipelines to natural gas pipelines, or the reinforcement 
of specific pipeline sections. 

In this paper, we address the Gas Transmission Network Expansion Planning (GTNEP) 
problem where the goal is to fulfill projected future gas contracts and to increase the 
reliability of a gas system under steady-state conditions. A Mixed-Integer NonLinear Pro¬ 
gramming (MINLP) formulation is proposed to model the design requirements and mini¬ 
mize expansion costs. Given the non-convex nature of the problem, a convex mixed-integer 
second-order cone relaxation is introduced. The proposed convex relaxation is based on 
four key ideas: (I) the introduction of variables for modeling the flux directions; (2) exact 
McGormick relaxations; (3) on/off constraints; and (4) valid integer cuts. Experimental 
results on the Belgian gas network and a test bed of large-scale synthetic instances demon¬ 
strate three key findings: 
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1. The convex relaxation produces tight lower bounds with high computational effi¬ 
ciency; 

2. The solution to the convex relxation can almost always be used to derive high-quality 
solutions to the original problem, leading to provably tight optimality gaps and, in some 
cases, global optimal solutions. 

3. The proposed approach scales to large-scale instances. 

The rest of this paper is organized as follows. Section presents the literature review. 
Section introduces the problem formulation. Section specifies the convex relaxation. 
Section presents the computational results and Section concludes the paper. 


2. Literature review 

The last four decades have seen an interest in natural gas planning problems such as 
optimal design, optimal reinforcement, and optimal expansion of gas pipeline systems. 
Algorithms for these problems can be classified in a number of different ways such as exact 


approaches ( 

Andre et al. 

2009, Bonnans et al. 

2011, Edgar et al. 

1978 

2001b, Wolf 

2004) 

and heuristics 

Andre 2010, Boyd et al 

1994 

Humpola et al. i 

'015 

Andre et al. 

2009 

Humpola and Fiigenschul 

^1 2014a) 

. Exact methods include cutting 

plan 

) and 

3S (Atamturk||2002 

Humpola and Fiigenschuhl 2014b 

Humpola et a 

1. 2015a, Poss|2011 

branch-and-bound 

) they use a variety 

(Andre||2010 

Elshiekh et al.||2013 

, Humpola and Fugenschuh||2011 

)) anc 

of commercial ( 

Bakhouya and De Wolf 

1008, I 

llshiekh et al. 2011 

L Soliman and Murtagh 


1982) and open-source (Pfetsch et al.||2012 Uster and Dilaveroglu 2014) solvers. Like this 


paper, much of the literature relies on approximations and relaxations to improve the 
tractablity of the underlying planning problems. Examples include continuous relaxations 


of the discrete design variables ( 

De Wolf and Smeers 

1996 

Hansen et al. 

1991 

Soliman 

and Murtagh 

1982 

) and approximation or relaxations of constraints ( 

Babonneau et al. 


Bakhouya and De Wolf||2008[ Humpola and Fugenschuh||2015 

Poss||2011 

). Common 


approaches for implementing these approximations/relaxations include succesive linear 


programming (De Wolf et al. 

1991, Hansen et al. 1991, O’Neill et al. 1979, Wilson et al. 

1988) and piecewise linearizations (Correa-Posada and Sanchez-Martm 2014, Markowitz 

and Manne||1957, Vajda 1964 

Zheng et al. 2010). 
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The contribution of this paper is a novel Second-Order Cone (SOC) relaxation that 
efficiently addresses the design of large-scale cyclic networks for which flow directions are 
unknown. The model captures physical, operational, contractual, and on/ofF constraints 
and includes models of regular pipelines, valves, short pipes, control valves, compressor 
stations, and regulators. Its dual solutions can almost always be converted to high-quality 
or optimal primal solutions. To the best of our knowledge, this combination of features 
has not appeared in the literature. Our paper focuses on the cost of building the network 
but can be generalized to include operational costs as well. 

We now provide an in-depth review of the most relevant works in the area of natu¬ 
ral gas expansion planning problems. One of the earliest papers that addresses natural 


gas design problems is (Edgar et al. 1978). It focuses on the optimal design of gunbar- 
rel and tree-shaped networks. Their objective minimizes the yearly cumulative operational 
and investment costs. The optimization variables include pipeline diameters, compression 


ratios, and the number of compressors. In their later work, Edgar et al. (2001b) present a 
MINLP formulation for the optimal design of a gas transmission network where the num¬ 
ber of compressor stations, the length and diameter of the pipeline sections, and the inlet 
and outlet pressures at each stations are optimized. They solve a simplified version of the 


problem in GAMS (GAMS Development Gorporation 2008) for a small instance (Edgar 
et al.||2001a ). 


Hansen et al. (1991) and Soliman and Murtagh (1982) propose a continuous relaxation 


for the network design problem. While Hansen et al. (1991) apply a successive linear pro¬ 
gramming method where a linear subproblem is solved to adjust the discrete choice of diam¬ 


eters, Soliman and Murtagh (1982) apply the commercial NLP solver MINOS Murtagh 


and Saunders (1998) to handle the relaxed subproblem. O’Neill et al. (1979) and Wil¬ 


son et al. (1988) focus on a problem where integer variables are used for the operational 


state of compressor stations and they also implement a method based on successive linear 
programming to solve the problem. 


De Wolf and Smeers (1996) address the optimal dimensioning of a known pipe network 


topology with an objective that combines the cost of purchasing gas and the capital expen¬ 
ditures for expansion. The authors formulate the problem as a continuous NLP that selects 
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pipeline diameters and solves the problem by means of a local optimizer. Based on this 


problem, Wolf (2004) derives conditions under which this problem is convex. Through the 
use of variational inequality theory, they show convexity of the nonlinear gas flow system 
under the assumption that the gas net inlet (pressure) is hxed at all supply and demand 


nodes. Bakhouya and De Wolf (2008) also present a case study on the same problem 
with separable transportation and gas objectives that leads to a two-stage problem for¬ 
mulation. In addition to design variables for the optimal pipe diameters, the authors add 
investment variables representing the maximal power of compressor stations to balance 
the pipeline construction costs and capital expenditures for increasing power in the com¬ 
pressor units. The authors hnd an initial solution by solving a convex problem where all 
pressure constraints are relaxed. Then, the complete problem is locally solved by means of 
the GAMS/CONOPT solver. In these works, numerical experiments are primarily focused 
on the Belgian gas transmission network. 


Andre et al. (2009) present a MINLP model to solve the investment cost minimization 


problem for an existing gas system that includes pipelines and regulators and omits com¬ 
pressor stations. The goal is to identify a set of pipeline sections to reinforce and to select 
an optimal diameter size for these sections based on a discrete set of diameters. Under 
the assumption that the network is radial (the head loss equations are convex when flows 
are hxed), the authors propose a continuous relaxation of the pipe diameters (continuous 
intervals). A branch-and-bound approach for a unique maximal demand scenario is applied 
to a segment of the French high-pressure natural gas transmission system. A complete 


review and extensions of these hndings are provided in (Andre 2010). 


Babonneau et al. (2009 2012) focus on the design and operation of a natural gas trans¬ 


mission system while minimizing investment, purchase, and transportation costs. The 
authors propose an approach based on a minimum energy principle that transforms the 
non-linear non-convex optimization problem into a convex problem. The underlying con¬ 
vex, bi-objective formulation is an approximation of the investment cost function and the 
cost of energy to transport the gas. Their continuous formulation is applied to non-divisible 
constraints such as a limited number of available commercial pipe dimensions. 
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Bonnans et al. (2011) presents several problems that include the minimization of com¬ 


pressor ratios and the sum of operations and investment costs. The authors propose a 
global optimization technique that is based on the combination of interval analysis with 
constraint propagation. 


Zheng et al. (2010) discusses different optimization models in the natural gas industry, 


including the compressor station allocation problem, the least gas purchase problem and 
optimal dimensioning of gas pipelines. The authors review solution techniques to solve 
the underlying models which include a piecewise linear programming algorithm and a 
branch-and-bound algorithm. 


Elshiekh et al. (2013) presents a model to optimize the design and operation of the 


Egyptian gas system, where continuous design variables for the length and diameter of 


pipelines are considered along with a modified Panhandle equation (Coelho and Pinho 


2007|. The complete model is directly solved by means of the computer-aided optimization 


software LINGO ( LINDO Systems]|1997 ). 

[Uster and Dilaveroglu (2014) address the cost minimization problem of designing a 
new natural gas transmission system and expanding an existing gas system. The authors 
propose a mathematical formulation to tackle the design/expansion network problem for 
a given multi-period planning horizon. The underlying MINLP model is formulated in 


AMPL and solved approximately with Bonmin (Bonami and Lee 2013). 


Humpola and Fiigenschuh (2014b) and Humpola et al. (2015a) present valid inequalities 


for a MINLP model of a design problem in gas transmission systems. Different relaxations 
are applied to the subproblems created after branching on the additive and design variables 
for the active and passive components. The resulting passive transmission subproblems, 
which are referred to as leaf problems, admit slack variables to independently relax the 
pressure domains and the flow conservation constraints. The proposed cutting planes aim 
at reducing the CPU time of a branch-and-cut-based outer approximation applied to the 


full model where construction costs are defined by a global constant. Atamturk (2002) 


and Poss (2011) also propose valid inequalities to reinforce the relaxation approach to the 


network design structure. 
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Humpola and Fiigenschuh (2015) examines different (convex) relaxations for snbprob- 


lems created while applying a branch-and-bound technique to a nonlinear network design 
problem. Cutting planes on the nonlinear potential loss constraints are used to strengthen 
the relaxed subproblems. 

Pfetsch et al.| (2012) focuses on the validation of nomination problem while considering 


regular pipes and valves, control valves, compressors and regulators. The authors describe 
a two-stage approach to solve the resulting MINLP problem and propose several modeling 
techniques and approaches to account for, e.g., pressure losses. They also developed several 


large test cases (GasLib 2014). These problems form the basis for many of the problems 
we consider in this paper. 


3. Problem Formulation 

This section derives the problem formulation (as a disjunctive program) in stepwise refine¬ 
ments. It starts by deriving a disjunctive formulation that is then refined by introducing 
flux variables. 


3.1. The Disjunctive Formultion 

Gas dynamics along a pipe is described by a set of partial differential equation (PDE) 


with both spatial and temporal dimensions (Osiadacz 1987, Thorley and They 1987, Sar- 


danashvili 2005): 


dtp + djpv) = 0, 

-i^pv\v\-gsuiOLp, 

(1) 

dt{pv) + djpv'^) + d^p = 

(2) 

p = ZRTp. 


(3) 


Gas velocity v, pressure p, and density p are defined for every point x along the pipe and 
evolve over time t. Z represents the gas compressibility factor, T the temperature, and R 
the gas constant. 

Equation ([^ enforces mass conservation, Equation ([^ describes momentum balance, 
and Equation ([^ defines the ideal gas thermodynamic relation. In Equation ([^, the first 
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term on the right-hand side represents the friction losses in a pipe of diameter D with 
friction factor f. The second term accounts for the gain or loss of momentum due to 
gravity g if the pipe is tilted by an angle a. In practice, frictional losses dominate the 
gravitational term which is dropped. One can also safely assume that the temperature does 
not fluctuate signihcantly along a pipe. If temperature gradients are significant, a spatial 
decomposition, splitting the pipe into temperature stable segments, can be adopted. 

Taking into account these assumptions, Equations Q,Q, and Q are rewritten in terms 
of pressure p and mass flux ({) = pv. 


dtp = -ZRTda,(l), 

(4) 

dxP = m\, 

(5) 


In this work, we assume that the system has reached a steady state after its first com¬ 
missioning and hence all time derivatives are set to zero. Given this assumption, a Graph 
Transmission Network (GTN) is represented by a graph Q = {J\f ,A) where M denotes the 
set of nodes representing connection points and A denotes the set of arcs. An arc is a 
triplet {a,i,j) consisting of a unique identifier a linking nodes i and j. For convenience, 
such a triplet {a,i,j) will be denoted by in the following. Observe that parallel arcs 
can link the same pair of nodes, e.g., we have arcs and a*^ in a GTN where a and a* 
are the unique identifiers of these arcs. 

By setting the time derivatives to zero, the total gas mass flux along a pipe a^- becomes 
constant, i.e., = (j)a- Hence Equations Q and ([^ simplify to 

Pl-P^j=Wa(l)a\(l)a\- ( 6 ) 

where . 

Gas System Components The problem formulation considers pipes, compressors, short 
pipes, resistors, and valves. Gompressors, short pipes, and valves are modelled as lossless 
pipelines, i.e., Wa = 0. A compressor installed on arc a^- can increase/decrease the pressure 
ratio Oa =Pj/Pi, within the bounds and q:“, where = 1 and q:“ > 1 is typical for 
most compressors. A bi-directional compressor can perform compression based on the flux 
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direction, i.e., it is able to invert the ratio to =Pi/pj if the flux is going from j to i. A 
standard valve features a binary on/off switch and a control valve has a continuous switch¬ 
ing mechanism to adjust pressure. Thus a valve installed on arc can increase/decrease 
the pressure ratio =Pj/pi, within the bounds and q;“, where > 0 and a" < 1 is 
typical for most control valves and = q:“ = 1 for all valves. Finally, a resistor is modelled 
as a pipeline with a particular (small) loss coefficient {w). 

Expansion Variables The set of arcs A = Ae U An includes existing arcs Ae = 
Pe U Ce U Ve, as wcll as new arcs An = 'Pn U C„. In this notation, Ve denotes the 
set of installed pipelines, resistors, and short pipes. Cg and Ve denote the set of existing 
compressors and valves (control and regular) respectively. Vn and C„ denote the set of 
new pipelines and new compressors respectively. A binary variable is assigned for each 
new pipe in Vn to model the expansion decision, i.e., z^ = l if pipeline a^j is installed 
and z^ = 0 otherwise. Variables z^^ {{i,j) ^ C„) have an equivalent interpretation for new 
compressors. A binary variable Va is used to control the switching actions of valves. 

Disjunctive Formulation Since the pressure variables only appear in a square form, the 
formulation uses the variable substitution Pi =pf (i G Af). Equations can be written as 

Pi - Pj = Wa4)a\4>a\ (a*j G V ( 7 ) 


Figure illustrates the curve of the function f{x,y) = y — mx|x| defined by the pressure 
drop equation Q. 

Since bi-directional compressor constraints depend on the flux direction, they can only 
be modelled using on/off or disjunctive constraints ( jHijazi et al. 2010, 2012). i.e., 

f PiO-i < Pj < PiOLl, if (/>a > 0 

1 PjOl\ <Pi< PjOcl, if Pa < 0. 


( 8 ) 


Given a set of injection (resp. demand) nodes X (resp. with mass flux injec¬ 

tion/demand Qi, the problem consists in finding an assignment of the expansion variables 
zP {oij G Vn), node pressures Pi {i G A7), and edge flows pa {chj G A), satisfying the Wey¬ 
mouth equations Q, the compressor constraints Q, and the following node conservation 
constraints: 

pa= pa+Qtii^Af) 
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Figure 1 The Gas Flow Equation f{x,y) = y — 'wx\x\. 

where qi = 0 for all i e N\{I LI V). Note that, in the steady-state model, injections are 
balanced, i.e., objective is to minimize the cost of expansion: 

min ^ Cazl+ ^ 

where Cq represents the cost of installing a new pipeline. The disjunctive formulation 
of the problem incorporating these ideas is presented in Model where /3* = (ct-)^ and 

/3r = «r- 

3.2. The Formulation Based on Flux Direction Variables 

This section presents a second formulation using flux direction variables to account for 
the disjunctive nature of the constraints. For every arc G A, we introduce two binary 
variables and y~ €{0,1} with the following semantics: y'l = 1 (resp. y~ = 1) if the flux 
moves from i to j (resp. from j to i) and 0 otherwise. The mass flux direction is captured 
by the following system of constraints: 

(1 - 2/a ) E Qfc < < (1 - y”) E 

< (1 -y+)y< A -Pj < (1 - 

. ya+ya=l- 

The first constraint ensures that y^ = 1 (resp. y~ = 1) if and only if 0a > 0 (resp. 0a < 0). 

Note that ^ Qk is an upper bound to the mass flux in a pipe. The second constraint 
fcGX 
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Model 1 The Disjunctive Formulation of the GTNEP. 

variables: 

J £ [Pi , PT ] \/i^N - squared pressure level variables 
^ A - mass flux on pipe (i,j) 

G {0, 1} Voy £Vn - binary expansion variables for pipes 
z'^ G {0, 1 } ya^j G - binary expansion variables for compressors 
Va G {0, 1 } G CVe U Ve - binary switch variables for valves 

objective: 

min CaZ^a+ (9a) 

subject to: 

E E ^a + Qi, VzGAA (9b) 

^i '^a4^a 1 4^a \ ? ^ij 


ZalPi PP "^a^al^al; G 


(9d) 

PiCtl < Pj < if Pa > 0, ya^j G Ce, 


(9e) 

P]OLa < Pi < PjOcl, if Pa < 0, VOy G Ce, 


(9f) 

PiOLa < Pj < if 0a > 0 and zl = 1, 

\f CLj^j G 

(9g) 

PjOL^ <Pi< PjOil, if 0a < 0 and z^ = 1, 

\f Qij G 

(9h) 

0a = 0 if z)) = 0, VOy G C„ 


(9i) 

< Pj < PiCC, if 0a > 0 and = 1, 

(Xlj G 

(9j) 

PjOpa <Pi< PjOiP if 0a < 0 and < = 1, 

VG-jj G Ve, 

(9k) 

0a = 0 if Ua = 0, ya^J G Ce 


(91) 


enforces a similar condition for the pressure difference. Using the variables and constraints, 
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the pressure drop equation can now be written without absolute value as 

{yl-Va) iJ- Pj) = Wa(l)l 

and the bi-directional compressor constraints are written as 


A«i-(i-2/^ 


;)<A 

< Ao^a + (1 

-ytm 


\f Qjj^j G (^e 

( 10 ) 


)(/3>i-/3 

D<A 

< A'O^a + (1 

-y-aW 


Vftjj G (Jg 

( 11 ) 

- (2 - ya 


-p]) 

< A < A^a 

+ (2-y+ 

-A)(A“ 

-A^a), Va,, gC„ 

( 12 ) 


-A)(/3>! 

.-(Ji) 

< A < A ^a 

+ (2-ya- 

-A)(A" 

-f3]aa), Va,, gC„ 

(13) 

Va, 

G C 

ij ^ ^n 





(14) 


iei iei 

The bi-directional valve constraints are written as 


Aai - (2 - - /3') < /5, < Aa“ + (2 - Va,, G V, (15) 

- (2 - l/a - Va){f3"a[ - (3l) < A < /3,q:“ + (2 - y" - - /3'a“), Va,, G Ve (16) 

- fa ^ < 0a < fa G Ve (17) 

iex iex 


The complete Mixed-Integer NonLinear Programming (MINLP) formulation based on flux 
direction variables is summarized in Model [H The continuous relaxation of Model [2] is 


non-convex due to Constraints (24c)-(24d). 


4. A Convex Relaxation of the GPNEP 

This section introduces a new mixed-integer second-order cone relaxation for Model 


4.1. The Variables 

For every pipe a^- G V, the relaxation introduces the auxiliary variable 7 ^ representing the 


product in Equations (18c)-(18d), i.e., 


7a = {yt - ya) (A - A)> ^ 


(19) 
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Model 2 The MINLP Formulation of the GTNEP. 

variables: 

A e [f3l,PT] Mi^N - squared pressure level variables 
(/)a G M £ A - mass flux on pipe (i,j) 

G {0,1} Voy £V„ - binary expansion variables for pipes 
G {0,1}, Vajj G A - binary flux direction variables 
G {0,1} Vajj G - binary expansion variables for compressors 
Va G {0,1} G CVe U Ve - binary switch variables for valves 


objective: 

min CaZ^,+ Y, (18a) 

subject to: 

E E (j)a + qi,^i£M (18b) 

{yl - yl) (A - A) = Va^ g Ve (isc) 

{yt - vl) (A - A) = G Vn (18d) 

(18e) 

fcGX fcGX 

(1 - yt) A < A - A < (1 - y-a) AA Va., G V (18f) 

@-[IZl) (18g) 

ya +ya =1: VOy-GA (18h) 


This product is then linearlized by a standard relaxation introduced by McCormick (1976) 
for bilinear functions, i.e., 


7a>A-A + (A-A“) (ya-ya +1) 

(20) 

7a>A-A + (A“-/3)) (y+-ya"-i) 

(21) 

7a<A-A + (A“-/3)) (y+-ya" + i) 

(22) 
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(23) 


This linearization is exact, since (y^ ~ya) take only discrete values. 


4.2. The Constraints 


The non-convex constraints (18c) can now be relaxed into 


7a >'Wa(l)l (Oy GPe). 


The on/off constraints (18d) represent another challenge for convexifying Model These 
constraints can be written as 


7a = ■Wa(t>l if = 1 (Op- e Vn) 


with a disjunctive second-order cone relaxation defined as 


7a > Wa4)l, if ^^a = 1 («0 ^ '^n)- 


Perspective formulations introduced by Hijazi et al. (2012) can be used to formulate the 
convex hull of such on/off constraints, giving the following rotated second-order cone con¬ 
straint: 


zl'ya>Wa4>l, G . 


The complete Mixed-Integer Second-Order Cone Programming (MISOCP) relaxation is 
presented in Model 

4.3. The Integer Cuts 

The MINLP and MISOCP formulations presented in Models and can be strenghtened 
by introducing the following valid integer cuts: 


T.ya+Yl ya“>l, ViGX (25) 


2/a > Vz G P 


( 26 ) 
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Model 3 The MISOCP Relaxation for the GPNEP. 
variables: 

A e [f3l,PT] Mi^N - squared pressure level variables 
(/)a G M £ A - mass flux on pipe (i,j) 

G {0,1} Voy £V„ - binary expansion variables for pipes 
G {0,1}, Vajj G A - binary flux direction variables 
7 o G M'*’ £V - auxiliary variables for bilinear terms 

G {0,1} Vajj G - binary expansion variables for compressors 
Va G {0,1} G CVe U Ve - binary switch variables for valves 


objective: 

min CaZ^a+ (24a) 

Q^ij^Pn 

subject to: 

E E yi£j\f (24:b) 

7a >Wa(l)l, VOy G Ve (24c) 

zlla > G Vn (24d) 

VOy-GA (24e) 

iGl iGX 

(1 - yt) (3l < A - A < (1 - 2/a ) Pr, Va,, G V (24f) 

@-[T3) (24g) 

2/a +2/a =b Va„GA (24h) 


Constraints (25) are generated for each injection node i£l: They state that at least one 
connected arc has an outgoing flow, taking the orientation of the arc into account to select 


the proper variables (?/+ for arcs leaving i and for arcs coming to i). Constraints (26) 
follow the same reasoning for demand nodes i£T>. 
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For a node i with degree two and no injection/demand {qi 
cut is valid 

'yt = yl* if aj^,a*,,GA 
^ yt = ya* if 
ya=yl* if 
ja=ya* if 


0), the following integer 


(27) 


It can be easily derived using the flux conservation constraints (24b) stating that, for a 


node with degree two and zero injection/demand, the flux direction of the incoming arc 


determines the flux direction of the outgoing arc. 


Finally, we can derive integer cuts for parallel pipelines; 


yt* =2/a : yaij,a*^GA. 


(28) 


Equations (28) state that parallel pipelines share the same flow direction. The validity of 


this cut follows from the pressure drop equations (18c) and the fact that parallel pipelines 
share the same pair of pressure variables. 


4.4. Converting the Convex Relaxation in a Feasible Solution to the GPNEP 

The solution to the relaxed Model is not always feasible for Model To obtain a feasible 
solution, we fix all the binary variables and use a nonlinear optimization solver to find a 
(locally) optimal solution to the resulting problem. When the local solver does not converge 
to a feasible solution, we consider primal solutions obtained when solving Model and 
repeat the process. 


5. Computational Experiments 

This section studies the performance of the proposed MINLP and MISOCP models and 


compares them with a model using a piecewise linear approximation. Section 5.1 describes 


the benchmarks and Section 5.2 the experimental setting and the various algorithms used. 


Section |5.3| and 5.4 report the computational results on the Belgian network and larger 


networks respectively, while Section 5.5 reports on the importance of the integer cuts. 
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Ref 


Network configuration 


|AA| 


\V\ 

Base 

New 

\Ve\ 

|Ce| 

\'Pn\ 

|Cn| 

A 

20 

6 

9 

24 

3 

0 

0 

Al 

22 

6 

9 

24 

3 

4 

2 

^2 

25 

6 

9 

24 

3 

7 

4 

^3 

29 

6 

9 

24 

3 

12 

5 

Bi 

20 

6 

9 

0 

0 

135 

12 

B 2 

20 

6 

9 

0 

0 

135 

12 

B 3 

20 

6 

9 

0 

0 

135 

12 

B 4 

20 

6 

9 

0 

0 

135 

12 


Table 1 Test Instances Based on the Belgian Network. 


5.1. The Benchmarks 

5.1.1. The Belgian Network Table shows the list of test instances based on the 
Belgian network depicted in Figure The table shows, for each benchmark, the number of 
nodes, sources, terminals, base pipelines, and compressor stations, as well as the number of 
new components (pipelines and compressors) that can potentially be added to the network 
topology. Note that benchmark A in Table is the real Belgium gas transmission network 
and Table shows the node characteristics for this 20-node, 24-pipeline, 3-compressor 
network. The reader is referred to the appendix of ( De Wolf and Smeers||2000 ) for further 
details on this network. Instances Ai — A^ captures various possible expansions to this base 
network. Figure and Tables and depict the location of the potential expansion plans 
and their associated data. The network expansion plans were designed for the Belgian gas 
network in order to capture events such as increase of the number of nominations and 
forecasting demand at the city gates, as well as excessive stress of the available supplies at 
the sources. 

Instances Bi — B 4 are based on the “optimization from scratch” benchmarks from 
( De Wolf and Smeers||2000 ) and ( Baboimeau et alT]|2012 ) (a = 1,1.6,5,ond6, respectively). 
In these papers, the authors use the Belgian gas network for a variation of the GTNEP 
problem which considers integrated functions of the gas merchant and transportation pro¬ 
cess. These benchmarks specify minimum and maximum production levels (see Table |^. 
Since the GTNEP assumes known gas nomination and production profiles, we computed 
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Figure 2 The Belgian Gas Network Base Configuration (case A). 


load and compression profiles based on optimal pressures provided in ( jBabonneau et ^ 
2012). Our instances also employ the same cost function as in ( Babonneau et ah]|2012 ) to 
compute the associated costs for building new pipelines, i.e., 

(l.04081"®Ilf/ +11.2155) 


where Hy and are the diameter and length of pipeline (i,j) respectively. (Babonneau 


et al. 2012) assumed continuous diameter choices. However, we used a discrete diameter 


values corresponding to the solution of ( De Wolf and Smeers||2000 ) and Table 4 of ( [Babon- 
neau et al. 2012). For completeness, the diameter choices are described in Table Note 
that the exclusive-set constraint is slightly different for these cases due to the existence of 
pre-defined parallel pipes. Within in each row of Table the solution must contain one 
and only diameter choice, and each set of parallel pipes must choose diameters from the 
same column of Table [2j 
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Pipe 

D, 

D2 

Ds 

D4 

^5 

(1,2) A 

890.0 

650.3 

610.8 

524.7 

512.1 

(1,2) B 

890.0 

650.3 

610.8 

524.7 

512.1 

(2,3) A 

890.0 

834.7 

784.0 

673.5 

657.3 

(2,3) B 

890.0 

834.7 

784.0 

673.5 

657.3 

(3,4) 

890.0 

998.9 

938.3 

806.0 

786.7 

(5,6) 

590.1 

604.3 

567.6 

487.6 

475.9 

(6,7) 

590.1 

0 

X 

X 

X 

(7,4) 

590.1 

671.7 

630.9 

542.0 

529.0 

(4,14) 

890.0 

829.9 

779.5 

669.7 

653.6 

(8,9) A 

890.0 

902.8 

848.0 

728.4 

711.0 

(8,9) B 

395.5 

902.8 

848.0 

728.4 

711.0 

(9,10) A 

890.0 

902.8 

848.0 

728.4 

710.9 

(9,10) B 

395.5 

902.8 

848.0 

728.4 

711.0 

(10.11) A 

890.0 

787.6 

739.8 

635.5 

620.1 

(10.11) B 

395.5 

787.6 

739.8 

635.5 

620.4 

(11,12) 

890.0 

979.8 

920.3 

790.6 

771.6 

(12,13) 

890.0 

915.1 

859.6 

738.4 

720.7 

(13,14) 

890.0 

952.6 

894.7 

768.6 

750.1 

(14,15) 

890.0 

1201.0 

1128.0 

969.0 

945.8 

(15,16) 

890.0 

1038.4 

975.3 

837.9 

817.7 

(11,17) 

395.5 

469.0 

440.5 

378.4 

369.3 

(17,18) 

315.5 

469.0 

440.5 

378.4 

369.3 

(18,19) 

315.5 

469.0 

440.5 

378.4 

369.3 

(19,20) 

315.5 

448.9 

421.7 

362.2 

353.5 


Table 2 


Pipe diameter choices from Table 4 of (jBabonneau et al.||2012 i 


5.1.2. Larger Networks Table [^describes the main data points for the larger bench¬ 
marks. Instance D is a real-life network case whose data is restricted for confidentiality 
reasons and we are not allowed to disclose its map or load profile. Instances E,F and G 
are part of a German network whose data, including the network configuration, maps, and 
load profiles, can be found in (Pfetsch et al.|j20I2). 


5.2. The Algorithms and the Experimental Setting 

This section reports computational results for three approaches: 

1. The MINLP formulation of the GTNEP as shown in Model [2| 

2. The MISOGP relaxation of the GTNEP as shown in Modelfollowed by the conver¬ 
sion presented in Section 


4.4 
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Node (Loc.) 

TypePl L 

(Loads) 

L L 

(Pressure) 
P P 

1 (Zeebrugge) 

I 

8.87 

11.594 

10.911288 

0 

77 

2 (Dudzele) 

I 

0 

8.4 

8.4 

0 

77 

3 (Brugge) 

V 

—oo 

-3.918 

-3.918 

30 

80 

4 (Zomergem) 


0 

0 

0 

0 

80 

5 (Loenhout) 

I 

0 

4.8 

2.814712 

0 

77 

6 (Antwerp) 

V 

-oo 

-4.034 

-4.034 

30 

80 

7 (Ghent) 

V 

-oo 

-5.256 

-5.256 

30 

80 

8 (Voeren) 

I 

20.34 

22.01 

22.012 

50 

66.2 

9 (Berneau) 


0 

0 

0 

0 

66.21 

10 (Liege) 

V 

-oo 

-6.365 

-6.365 

30 

66.2 

11 (Warnand) 


0 

0 

0 

0 

66.2 

12 (Namur) 

V 

-oo 

-2.12 

-2.12 

0 

66.2 

13 (Anderlues) 

X 

0 

1.2 

1.2 

0 

66.2 

14 (Peronnes) 

X 

0 

0.96 

0.96 

0 

66.2 

15 (Mons) 

V 

-oo 

-6.848 

-6.848 

0 

66.2 

16 (Blaregnies) 

V 

-oo 

-15.616 

-15.616 

50 

66.2 

17 (Wanze) 


0 

0 

0 

0 

66.2 

18 (Sinsin) 


0 

0 

0 

0 

63 

19 (Arlon) 

V 

-oo 

-0.222 

-0.222 

0 

66.2 

20 (Petange) 

V 

-oo 

-1.919 

-1.919 

25 

66.2 


Table 3 The Belgian Gas Network Data from ( [Pe Wolf and Smeers||2000^ . This data is used for the 
A Problems, f- On Problems Al, A2, and A3,the pressure bounds are [0,59.851968], [0,59], and [0,59.85] 

respectively. 


3. A MIP formulation based on a Piecewise Linear Approximation (PLA-MIP) of the 


quadratic functions; The PLA-MIP formulation follows the derivation in (Correa-Posada 


and Sanchez-Martin 2014, De Wolf and Smeers||2000 ) and uses 60 segments. 


All the experiments weere conducted on a computer with two Intel Xeon CPU X5670 
processors (2.93GHz) with 6 cores each. The computer has 64 GB DIMM 1333MHz RAM 
and runs the Ubuntu 14.04 LTS operating system. The MINLP formulation is solved using 
SCIP 3.1.1 ( AchterbergI 2009) compiled with Ipopt 3.12.3 and Cplex 12.6. The PLA-MIP 


formulation is solved using CPLEX 12.6 (ILOG CPLEX Optimization Studio 2013). The 
MISCOP formulation is solved with CPLEX 12.6 and the conversion is performed by 


IPOPT 3.12.3 (Wchter and Biegler 


2006). 
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Load (L) profiles (MMscf) 


Node 

B, 

B 2 

Bs 

Bi 

1 

9.5883 

9.8225 

9.8218 

9.7205 

2 

8.1833 

8.3447 

8.1340 

8.3628 

3 

-3.9180 

-3.9180 

-3.9180 

-3.9180 

4 

0.0000 

0.0000 

0.0000 

0.0000 

5 

4.0315 

4.0432 

4.0383 

4.0364 

6 

-4.0315 

-4.0432 

-4.0383 

-4.0364 

7 

-5.2413 

-5.2644 

-5.2562 

-5.2644 

8 

22.012 

22.0120 

22.0120 

22.0120 

9 

0.0000 

0.0000 

0.0000 

0.0000 

10 

-6.4744 

-6.4951 

-6.3970 

-6.3816 

11 

0.0000 

0.0000 

0.0000 

0.0000 

12 

-2.1929 

-2.1191 

-2.1162 

-2.0984 

13 

1.2162 

1.3225 

1.0802 

1.1591 

14 

0.9840 

0.6164 

1.0776 

1.0235 

15 

-6.4056 

-6.5885 

-6.8366 

-6.8857 

16 

-15.6119 

-15.5904 

-15.4616 

-15.5899 

17 

0.0000 

0.0000 

0.0000 

0.0000 

18 

0.0000 

0.0000 

0.0000 

0.0000 

19 

-0.2059 

-0.2312 

-0.2269 

-0.2164 

20 

-1.9337 

-1.9112 

-1.9131 

-1.9236 


Table 4 The Load Profiles Computed from Optimal Pressures Provided in ( [Babonneau et aL]|2012^ . 

All compression ratios were derived as 1.0. 



Figure 3 The Belgian Gas Network Expansion Plans (Instances A 1 -A 3 ). 
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Node 

Town 

Lat. 

Long. 

p(*) 

p(*) 

21 

Instance Ai 
Bois 

50.400676 

5.855991 

14 

66 

22 

Koninklijke 

50.806672 

4.481877 

14 

66 

21 

Instance A 2 
Heist 

51.095651 

4.744616 

20 

70 

22 

Zoutleeuw 

50.858734 

5.115404 

20 

70 

23 

Beaufays 

50.552195 

5.670182 

20 

70 

24 

Gouvy 

50.231757 

5.966813 

20 

70 

25 

Ettelbruck 

49.861370 

6.073930 

20 

70 

21 

Instance A 3 
Jabbeke 

51.204699 

3.086440 

14 

66 

22 

Torhout 

51.072867 

3.118026 

14 

66 

23 

Kortrijk 

50.790711 

3.230636 

14 

66 

24 

Bois-de-Barry 

50.580151 

3.521773 

14 

66 

25 

Lobbes 

50.353208 

4.263261 

20 

70 

26 

Senzeille 

50.124840 

4.433550 

20 

70 

27 

Gedinne 

49.980230 

4.851030 

20 

70 

28 

Ghiny 

49.806832 

5.274004 

20 

70 

29 

Pigneule 

49.735878 

5.471758 

20 

70 


Table 5 


Locations of Nodes of the Expansion Plans for the Belgian Gas Network. These nodes do not 

have injections. 


5.3. Results on the Belgian Network 

Table shows the sizes of the underlying models in terms of the number of binary and 
continuous variables and the number of linear and quadratic constraints for each instance. 
Table presents the computational results and reports the CPU time in seconds and 
the upgrade cost in $ x 10^ for each approach. The computational results show that the 
MISOCP approach outperforms both the MINLP and the PLA-MIP and that the solution 
to the MISOCP always converts to a feasible and optimal solution. The PLA-MIP approach 
has both computational and accuracy issues, as it significantly underestimates the optimal 
objective value and is rather slow. 


The results for problems B1-B4 are interesting as the expansion costs are considerably 


lower than reported by 

Babonneau et al. 

(20U 

!) for the same operating conditions. In 

Table 4 of their paper. 

Babonneau et al. I 

2012 ) 

report expansion costs of 15669, 14252, 


11610, and 11274 for B1-B4. Their solutions are feasible and have the same operating 
cost as our model. Of course, it is important to note than their solutions were obtained 
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Node 

Node 

w 

c 

Instance Ai 



9 

21 

0.929 

67.19 

21 

18 

0.808 

77.26 

6 

22 

0.785 

79.50 

22 

14 

0.766 

81.44 

Instance A 2 



5 

21 

1.052 

59.29 

21 

21* 

Compressor 

1500.1 

22 

11 

0.967 

64.52 

8 

23 

1.933 

32.28 

23 

24 

0.876 

71.18 

24 

24* 

Compressor 

1500.1 

25 

19 

1.339 

46.59 

21* 

22 

0.980 

63.65 

24* 

25 

0.866 

72.08 

Instance A 3 



1 

21 

2.257 

27.65 

2 

21 

4.546 

13.73 

21 

21* 

Compressor 

1500.1 

22 

23 

1.121 

55.66 

23 

23* 

Compressor 

1500.1 

24 

15 

1.073 

58.14 

15 

25 

1.483 

42.09 

25 

26 

1.289 

48.40 

26 

26* 

Compressor 

1500.1 

27 

28 

1.010 

61.79 

28 

29 

2.232 

27.96 

29 

19 

1.423 

42.09 

21* 

22 

2.448 

25.50 

23* 

24 

1.165 

53.56 

26* 

27 

1.071 

58.28 


Table 6 Locations of Pipes of the Expansion Plans for the Belgian Gas Network. * denotes 
introduced dummy node for 0 length compressor arcs. 


through a model that minimizes operating and expansion costs, which conld make it harder 
to determine the best design for particular operating conditions. Still, this comparison 
highlights the strengths of the formulation proposed in this paper. 

5.4. Scalability Results 

We now stndy whether the resnlts on the Belgian networks continue to hold on larger 
instances. To assess scalability and robustness, we stress the networks by gradually increas- 
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Network configuration 


Base New 


Ref. 

|A| 


\V\ 

\'Pe\ 

|Ce| 

l^nl 

D 

60 

2 

24 

55 

4 

55 

E 

40 

3 

29 

39 

6 

39 

F 

135 

6 

99 

141 

29 

141 

G 

582 

31 

129 

609 

5 

278 


Table 7 Larger Instances of Gas Networks. 


MINLP PLA-MIP MISOCP 


Bench. 

BV 

CV 

LC 

QC 

BV 

CV 

LC 

QC 

BV 

CV 

LC 

QC 

A 

54 

49 

254 

96 

1494 

1837 

3931 

0 

54 

73 

398 

24 

Al 

70 

59 

320 

112 

1750 

2151 

4605 

0 

66 

91 

488 

28 

^2 

85 

69 

389 

124 

1945 

2391 

5120 

0 

78 

107 

575 

31 

^3 

103 

80 

463 

144 

2263 

2776 

5954 

0 

91 

128 

679 

36 

- 61 , 2 , 3,4 

354 

1154 

464 

357 

7314 

8737 

19067 

0 

238 

373 

1850 

116 


Table 8 The Sizes of the Mathematical models for Belgian Network Instances. (BV: Binary variables, 
CV: Continuous variables, LC: Linear constraints, QC: Quadratic constraints). 


Bench. 

MINLP 

PLA-MIP 

MISOCP 

CPU 

Obj 

CPU 

Obj 

CPU 

Obj 

A 

0.02 

0.0 

0.6 

0.0 

0.03 

0.0 

Al 

0.06 

144 

0.7 

144 

0.05 

144 

A 2 

0.06 

1687 

1.4 

187 

0.1 

1687 

As 

0.06 

1780 

1.9 

280 

0.06 

1780 

Bi 

1.89 

11181 

1089 

10353 

0.3 

11181 

B 2 

3.17 

11181 

1781 

10361 

0.6 

11181 

63 

3.53 

11181 

1538 

10352 

0.6 

11181 

64 

3.82 

11181 

1570 

10352 

0.3 

11181 


Table 9 Computational Results on the Belgian Network Instances: The Objective Value is in $ and 

the CPU Time in Seconds. 


ing the production and consumption levels from 5% up to 300% while considering solely 
the addition of a parallel pipe for each existing pipeline in the base configuration of the 
gas systems (i.e., |C„| =0). Table presents the sizes of the mathematical models. In all 
of these results, we denote whether or not the MISOCP and MINLP solutions are exact, 
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lower bounds, or upper bounds on the MINLP solutions. Lower bounds for the MINLP 
are also derived by subtracting the optimality gap from any primal feasible solution. 

Table [TT| presents the computational results on instance D which is based on proprietary 
natural gas network in the United States. Observe that the PLA-MIP model systematically 
underestimates the objective function and returns infeasible solutions. As we will see, this 
is systematic on all larger benchmarks. The MISOP approach returns optimal solutions 
for all but one case. Both the MINLP and MISOCP prove infeasibility of the most stressed 
network. 


Table 12 presents the computational results on instance E which is based on gaslib-40 
( GasLib||2014 ). The MISOP approach returns optimal solutions, or proves infeasibilities in 
all cases. The MISOP model is one order of magnitude faster than the MINLP model. 


Table 13 presents the computational results on instance F, which is based on gaslib- 


135 (GasLib 2014) and is particularly challenging. The MINLP approach finds optimal 
solutions up to the 25% case and spends considerable time doing so. It finds an upper 
bound to the 50% case but does not return any information on the 75% and 100% cases. 
In contrast, the MISOCP approach finds optimal solutions to the 0%, 5%, 25%, and 50% 
cases, all below 10 seconds. It finds lower bounds on the 75% and 100% cases reasonably 
fast. Both the MINLP and the MISOCP prove infeasibility of the three most stressed 
instances. 


Table 14 presents very interesting results for instance G, which is based on gaslib-582 


(GasLib 2014). The MINLP approach cannot find feasible solutions on any of the cases 
but the 300% case which is shown infeasible. Both the MINLP and PLA-MIP approaches 
have numerical issues with these problems. The MISOP approach finds optimal solutions 
up to the 50% case and for the 150% case and proves infeasibilities for the 200% and 300% 
cases. For the 75%-125% cases, the MISCOP times out but returns upper bounds to the 
optimal solution with duality gaps ranging from 7.65% to 51.3%. 

Overall, these results demonstrate the benefits of the MISOCP approach. The MISOCP 
approach almost always finds optimal solutions much faster than the MINLP when both 
return optimal solutions. It also finds optimal solutions or proves infeasibility in many case 
for the larger benchmarks, while the MINLP approach does not return feasible solutions. 
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MINLP PLA-MIP MISOCP) 


Ref. 

BV 

CV 

LC 

QC 

BV 

CV 

LC 

QC 

BV 

CV 

LC 

QC 

D 

283 

174 

1093 

440 

6883 

8330 

18018 

0 

228 

339 

1753 

no 

E 

207 

124 

792 

312 

4887 

5920 

12796 

0 

168 

241 

1260 

78 

F 

763 

446 

2886 

1128 

17683 

21430 

46304 

0 

622 

869 

4578 

282 

G 

2101 

1469 

8058 

2256 

35941 

44433 

92848 

0 

1823 

2311 

11442 

564 


Table 10 Size of the Mathematical Models: BV: Binary variables, CV: Continuous variables, LC: 

Linear constraints, QC: Quadratic constraints. 


Stresss 

level 

MINLP 

PLA-MIP 

MISOCP 

CPU 

Obj 

CPU 

Obj 

CPU 

Obj 

0 ^ 

0.1 

0 . 00 * 

3.0 

0.00 

0.1 

0 . 00 * 

5 % 

0.5 

3 . 50 * 

1.8 

0.00 

0.6 

3 . 50 * 

10 % 

1.6 

23 . 83 * 

12.2 

23.22 

0.5 

23 . 83 * 

25 % 

2.1 

92 . 24 * 

14.0 

83.99 

0.6 

92 . 24 * 

50 % 

1.5 

145 . 58 * 

14.8 

136.2 

0.5 

145 . 58 * 

75 % 

0.6 

191 . 80 * 

11.0 

184.0 

0.6 

191 . 8 * 

100 % 

3.0 

287 . 00 * 

12.5 

209.03 

0.7 

28 L 99 A 

125 % 

0.2 

t 

1.6 

t 

0.2 

t 


Table 11 Computational Results on Instance D. Obj: $ x 10^, CPU time: in seconds. Solution status: 
* = Proven optimal; A = Lower bound; v = Upper bound; f = Infeasible; J = Unknown. 


Stress 

level 


MINLP 


PLA-MIP 


MISOCP 


CPU 

Obj 

Gap 

CPU 

Obj 

Gap 

CPU 

Obj 

Gap 

0 ^ 

1.6 

0 . 00 * 

0.0 

10.2 

0.00 

0.0 

0.2 

0 . 00 * 

0.0 

5 % 

6.3 

11 . 92 * 

0.0 

23.5 

0.00 

0.0 

0.7 

11 . 92 * 

0.0 

10 % 

6.8 

32 . 83 * 

0.0 

20.6 

0.00 

0.0 

0.4 

32 . 83 * 

0.0 

25 % 

5.6 

41 . 08 * 

0.0 

30.9 

32.8 

0.0 

0.6 

41 . 08 * 

0.0 

50 % 

8.1 

156 . 06 * 

0.0 

11.5 

32.8 

0.0 

0.9 

156 . 06 * 

0.0 

75 % 

12.0 

333 . 01 * 

0.0 

21.8 

121.1 

0.0 

0.7 

333 . 00 * 

0.0 

100 % 

12.1 

551 . 64 * 

0.0 

17.5 

122.37 

0.0 

0.8 

551 . 64 * 

0.0 

125 % 

2.2 

t 

- 

33.0 

256.22 

0.0 

0.4 

t 

- 

150 % 

0.8 

t 

- 

27.6 

t 

- 

0.3 

t 

- 


Table 12 Computational Results on Instance E. Obj: $ x 10^, CPU time: in seconds. Solution status: 
* = Proven optimal; A = Lower bound; v = Upper bound; f = Infeasible; J = Unknown. 
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Stress 

MINLP 


PLA-MIP 


MISOCP 


level 

CPU 

Obj 

Gap 

CPU 

Obj 

Gap 

CPU 

Obj 

Gap 

0^ 

0.85 

olo^ 

0.0 

136.3 

0.0 

0.0 

1.3 

h0*“ 

0.0 

5% 

101.8 

0.0* 

0.0 

120.0 

0.0 

0.0 

1.0 

0.0* 

0.0 

10% 

36707.3 

15.04* 

0.0 

125.8 

0.0 

0.0 

2.4 

O.OA 

0.0 

25% 

457.9 

60.4* 

0.0 

124.4 

0.0 

0.0 

4.4 

60.4* 

0.0 

50% 

86962.9 

182.7V 

91.7 

166.7 

60.4 

0.0 

7.6 

95.3* 

0.0 

75% 

86933.9 


- 

119.8 

60.4 

0.0 

40.5 

451.5A 

0.0 

100% 

87334.2 


- 

119.5 

149.6 

0.0 

104.6 

1234.2A 

0.0 

125% 

6.8 

t 

- 

125.7 

149.6 

0.0 

1.8 

t 

0.0 

150% 

3.4 

t 

- 

206.7 

486.0 

0.0 

1.1 

t 

0.0 

200% 

0.4 

t 

- 

11.6 

t 

- 

0.4 

t 

0.0 


Table 13 Computational Results on Instance F. Obj: $ x 10^, CPU time: in seconds, Solution status: 
* = Proven optimal; A = Lower bound; v = Upper bound; f = Infeasible; J = Unknown. 


Stress 

MINLP 


PLA-MIP 


MISOCP 


level 

CPU 

Obj 

Gap 

CPU 

Obj 

Gap 

CPU 

Obj 

Gap 

0^ 

86400.0 

“T 

- 

62012.9 

6.87 

0.0 

2.7 

0.00* 

0.0 

5% 

86400.0 


- 

29655.1 

2.78 

0.0 

4.4 

0.00* 

0.0 

10% 

86400.0 


- 

86400.0 

4.65 

40.22 

21.1 

0.00* 

0.0 

25% 

86400.0 


- 

2153.2 

8.65 

0.0 

40.9 

0.00* 

0.0 

50% 

86400.0 


- 

3670.2 

t 

- 

164.0 

14.93* 

0.0 

75% 

86400.0 


- 

0.21 

t 

- 

86402.1 

111.99V 

51.3 

100% 

86400.0 


- 

5.31 

t 

- 

86401.6 

332.53V 

7.65 

125% 

86400.0 


- 

5.31 

t 

- 

86402.4 

524.82V 

11.74 

150% 

86400.0 


- 

5.29 

t 

- 

53321.3 

590.84* 

0.0 

200% 

86400.0 


- 

5.02 

t 

- 

16.7 

t 

- 

300% 

4.4 

t 

- 

0.12 

t 

- 

0.9 

t 

- 


Table 14 Computational Results on Instance G. Obj; $ x 10^, CPU time; in seconds, Solution status: 
* = Proven optimal; A = Lower bound; v = Upper bound; f = Infeasible; = Unknown. 


5.5. The Importance of Integer Cuts 


Table 15 describes the performance of the MISCOP on instances E, F, and G when the 
integer cuts are not used. As can be seen, the integer cuts, which were used both in the 
MINLP and MISOCP models, are critical to obtain an efficient MISOCP implementation. 


6. Concluding Remarks 

This paper considered the expansion of natural gas networks, a critical process involving 
substantial capital expenditures with complex decision-support requirements. It proposed 
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Stress 

level 

Instance E 


Instance F 


Instance G 


CPU 

Obj 

Gap 

CPU 

Obj 

Gap 

CPU 

Obj 

Gap 

0^ 

0.9 

0.00* 

0.0 

3310.9 

0.00* 

0.0 

242.2 

0.00* 

0.0 

5% 

1.8 

11.92* 

0.0 

83.5 

0.00* 

0.0 

14.2 

0.00* 

0.0 

10% 

2.7 

32.83* 

0.0 

120.7 

O.OOA 

0.0 

301.5 

0.00* 

0.0 

25% 

3.2 

41.08* 

0.0 

86419.5 

60.44V 

75.1 

86400.3 


- 

50% 

8.5 

156.06* 

0.0 

17693.1 

95.32^ 

0.0 

8271.05 

14.93* 

0.0 

75% 

6.7 

333.01* 

0.0 

86409.9 

451.59 

59.2 

86404.4 

111.99V 

79.4 

100% 

3.8 

551.64* 

0.0 

86404.5 

1234.23 

44.2 

87193.1 

332.32V 

87.9 

125% 

1.8 

t 

- 

90.9 

t 

- 

86401.8 

524.82V 

16.1 

150% 

1.0 

t 

- 

7.5 

t 

- 

86408.9 

245.80A 

58.5 

200% 

0.7 

t 

- 

2.0 

t 

- 

13318.6 

t 

- 

300% 

0.0 

t 

- 

0.0 

t 

- 

3.0 

t 

- 


Table 15 Computational Results on Instances E, F, and G without the Integer Cuts. Obj: $ x 10^, 


CPU time: in seconds, Solution status: * = Proven optimal; A = Lower bound; v = Upper bound; f = 

Infeasible; = Unknown 


a convex mixed-integer second-order cone relaxation for the gas expansion planning prob¬ 
lem under steady-state conditions in order to address the fact that state-of-the-art global 
optimisation solvers are unable to scale up to real-world size instances. The resulting MIS- 
OCP model offers tight lower bounds with high computational efficiency. In addition, the 
optimal solution of the relaxation can often be used to derive high-quality solutions to the 
original problem, leading to provably tight optimality gaps and, in some cases, global opti¬ 
mal solutions. The convex relaxation is based on a few key ideas, including the introduction 
of flux direction variables, exact McCormick relaxations, on/off constraints, and integer 
cuts. Numerical experiments are conducted on the traditional Belgian gas network, as well 
as other real larger networks. The computational results demonstrate that the MISOCP 
model is faster than the originating MINLP model by one or two orders of magnitude 
on the Belgian network instances. They also show that the MISOCP model scales well to 
large and stressed instances. 
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